**************************************************************
* ELECTRICITY SPILOVERS PROJECT
*  HOUSEHOLD SIMULATION MODEL
**************************************************************
set more off
frame create simulation
frame change simulation
local sim 10000 //Number of simulations

********************************************************************
********************************************************************
********************************************************************
*CLEANING RASS SURVEY DATA
********************************************************************
********************************************************************
********************************************************************
import delimited $data_main\Survdata.csv, clear

keep numi ctevpage shwrday bathsday showerhd aerators lndryeqp ///
	 cwhwld cwwwld cwcwld drylds dishwash dwloads sptyp sphtf pltyp ///
	 smrnflthr smdflthr snitflthr wmrnflthr wmdflthr wnitflthr numwash ///
	 plhtf smrnflthr wlwtrpmp czt24 res whelec edry ///
	 new_wht_uec new_swp_uec new_pmp_uec new_wpm_uec new_dwh_uec ///
	 new_cws_uec new_edy_uec new_spa_uec new_sph_uec

****	 
*Keeping only SoCal Single-Family Homes
* Note: - Have 2,188 SFHs in climate zone 9 (where Burbank is located)
*       - Could expand to all of coastal SoCal and get 8,442 SFHs
keep if res==1 
//keep if czt24==9
keep if czt24==5|czt24==6|czt24==7|czt24==9|czt24==10|czt24==14|czt24==16

****	 
*Cleaning and Labeling Variables
label var numi "Average No. Ppl in House"
label var ctevpage "No of Central Evap Coolers"
	replace ctevpage=. if ctevpage>3 //Missing observations
	replace ctevpage=1 if (ctevpage>0) & !missing(ctevpage)
label var whelec "Electric water heater"
label var shwrday "Showers per day"
	replace shwrday=. if shwrday>10 //Missing observations
label var bathsday "Baths per day"
	replace bathsday=. if bathsday>10 //Missing observations
label var showerhd "Owns low flow showerhead"
	replace showerhd=0 if showerhd==3 //Doesn't own
	replace showerhd=1 if showerhd==2 //Owns for some showers
	replace showerhd=. if showerhd>3 //Missing observations
label var aerators "Owns low flow showerhead"
	replace aerators=. if aerators>3 //Missing observations
	replace aerators=0 if aerators==3 //Doesn't own
	replace aerators=1 if aerators==2 //Owns some
label var lndryeqp "Owns laundry equipment"
	replace lndryeqp=0 if lndryeqp>=2 //Doesn't own
label var cwhwld "No. Hot Water Loads per Year"
	replace cwhwld=. if cwhwld>10 //Missing observations
	replace cwhwld=cwhwld*52 //Scaling up to loads per year
label var cwwwld "No. Warm Water Loads per Year"
	replace cwwwld=. if cwwwld>10 //Missing observations
	replace cwwwld=cwwwld*52 //Scaling up to loads per year
label var cwcwld "No. Cold Water Loads per Year"
	replace cwcwld=. if cwcwld>10 //Missing observations
	replace cwcwld=cwcwld*52 //Scaling up to loads per year
gen cwlds=cwhwld+cwwwld+cwcwld
	label var cwcwld "Clothes Washer Loads per Year"	
label var edry "Owns Electric Dryer"
label var drylds "No. Dryer Loads per Year"
	replace drylds=. if drylds>10 //Missing observations
	replace drylds=drylds*52 //Scaling up to loads per year
label var dishwash "Owns Dishwasher"
	replace dishwash=0 if dishwash==2 //Doesn't own
label var dwloads "No. Dishwasher Loads per Year"
	replace dwloads=. if dwloads>10 //Missing observations
	replace dwloads=dwloads*52 //Scaling up to loads per year
label var sptyp "Owns other spa"
	replace sptyp=. if sptyp>3 //Missing observations
	replace sptyp=0 if sptyp>=2  & !missing(sptyp) //Doesn't own
label var sphtf "Owns electric spa"
	replace sphtf=. if sphtf>6 //Missing observations
	replace sphtf=0 if sphtf>=2 & !missing(sphtf)  
	replace sphtf=0 if sptyp==0
	replace sptyp=0 if sphtf==1 //Setting spa ownership to zero if electric spa
label var pltyp "Owns pool"
	replace pltyp=. if pltyp>3 //Missing observations
	replace pltyp=0 if pltyp>=2  & !missing(sptyp) //Doesn't own
label var plhtf "Electric pool heater"
	replace plhtf=. if plhtf>7 //Missing observations
	replace plhtf=0 if (plhtf==1|plhtf==2|plhtf==5|plhtf==6|plhtf==7) & !missing(plhtf) //Other heat 	
	replace plhtf=1 if(plhtf==3|plhtf==4) & !missing(plhtf) //Electric heat 
label var wlwtrpmp "Owns electric water pump"
	replace wlwtrpmp=. if wlwtrpmp>2 //Missing observations
	replace wlwtrpmp=0 if wlwtrpmp==1 & !missing(wlwtrpmp)  
	replace wlwtrpmp=1 if wlwtrpmp==2 & !missing(wlwtrpmp) 
label var czt24 "Title 24 climate zone"
label var res "Building type"
label var numwash "No. Clothes Washes per week"
label var new_wht_uec "Water Heat Electric UEC"
label var new_swp_uec "Evap Cooler Electric UEC"
label var new_pmp_uec "Pool Pump UEC"
label var new_wpm_uec "Well Pump UEC"
label var new_dwh_uec "Dishwasher UEC"
label var new_cws_uec "Clothes washer UEC"
label var new_edy_uec "Electric dryer UEC"
label var new_spa_uec "Electric spa UEC"
label var new_sph_ue "Electric spa heat UEC"
 
*Pool filter use 
*Notes: - Want average pool filter use for year based on average daily use stats
replace smrnflthr=0 if smrnflthr==1 //No use
	replace smrnflthr=1.5 if smrnflthr==2 //1-2 hours use
	replace smrnflthr=3.5 if smrnflthr==3 //3-4 hours use
	replace smrnflthr=6 if smrnflthr==4 //5-7 hours use
	replace smrnflthr=10 if smrnflthr==5 //8-12 hours use
	replace smrnflthr=. if smrnflthr>5 //Missing observations
replace smdflthr=0 if smdflthr==1 //No use
	replace smdflthr=1.5 if smdflthr==2 //1-2 hours use
	replace smdflthr=3.5 if smdflthr==3 //3-4 hours use
	replace smdflthr=6 if smdflthr==4 //5-7 hours use
	replace smdflthr=10 if smdflthr==5 //8-12 hours use
	replace smdflthr=. if smdflthr>5 //Missing observations
replace snitflthr=0 if snitflthr==1 //No use
	replace snitflthr=1.5 if snitflthr==2 //1-2 hours use
	replace snitflthr=3.5 if snitflthr==3 //3-4 hours use
	replace snitflthr=6 if snitflthr==4 //5-7 hours use
	replace snitflthr=10 if snitflthr==5 //8-12 hours use
	replace snitflthr=. if snitflthr>5 //Missing observations	
egen slftuse=rowtotal(smrnflthr smdflthr snitflthr) 
	replace slftuse=slftuse*182.5 //Multiply by half year for 6 month use

replace wmrnflthr=0 if wmrnflthr==1 //No use
	replace wmrnflthr=1.5 if wmrnflthr==2 //1-2 hours use
	replace wmrnflthr=3.5 if wmrnflthr==3 //3-4 hours use
	replace wmrnflthr=6 if wmrnflthr==4 //5-7 hours use
	replace wmrnflthr=10 if wmrnflthr==5 //8-12 hours use
	replace wmrnflthr=. if wmrnflthr>5 //Missing observations
replace wmdflthr=0 if wmdflthr==1 //No use
	replace wmdflthr=1.5 if wmdflthr==2 //1-2 hours use
	replace wmdflthr=3.5 if wmdflthr==3 //3-4 hours use
	replace wmdflthr=6 if wmdflthr==4 //5-7 hours use
	replace wmdflthr=10 if wmdflthr==5 //8-12 hours use
	replace wmdflthr=. if wmdflthr>5 //Missing observations
replace wnitflthr=0 if wnitflthr==1 //No use
	replace wnitflthr=1.5 if wnitflthr==2 //1-2 hours use
	replace wnitflthr=3.5 if wnitflthr==3 //3-4 hours use
	replace wnitflthr=6 if wnitflthr==4 //5-7 hours use
	replace wnitflthr=10 if wnitflthr==5 //8-12 hours use
	replace wnitflthr=. if wnitflthr>5 //Missing observations	
egen wlftuse=rowtotal(wmrnflthr wmdflthr wnitflthr) 
	replace wlftuse=wlftuse*182.5 //Multiply by half year for 6 month use
	
egen plfltuse= rowtotal(slftuse wlftuse) 
drop smrnflthr smdflthr snitflthr wmrnflthr wmdflthr wnitflthr slftuse wlftuse
label var plfltuse "Yearly pool filter use (hours)"


********************************************************************
********************************************************************
********************************************************************
*ESTIMATING SATURATION RATES & UEC DISTRIBUTIONS
********************************************************************
********************************************************************
********************************************************************

*****
*Saturation rates
* Notes: - Based on average ownership of each electric-water appliance in sample
sum ctevpage //Own evaporative cooler
	local evap_sat=r(mean) 
sum whelec //Own electric water heater
	local wht_sat=r(mean) 
sum lndryeqp //Owns clothes washer
	local cws_sat=r(mean) 
sum edry //Owns electric dryer
	local edy_sat=r(mean)  	
sum dishwash //Owns dishwasher
	local dwh_sat=r(mean)  	
sum pltyp //Owns pool pump
	local pmp_sat=r(mean)  	
sum sptyp //Owns other spa
	local spa_sat=r(mean) 
sum sphtf //Owns other spa
	local sph_sat=r(mean)  
sum wlwtrpmp
	local wpm_sat=r(mean)

*****
*UEC Distributions
* Notes: - Assume all distributions are beta (ensures non-negative #s)
*        - To estimate: (i) Divide UEC by max value to ensure [0,1]
*                       (ii) Fit beta distribution, collect alpha and beta
*                       (iii) Draw from beta(alpha,beta) and multiply by max
*Evaporative cooler
sum new_swp_uec  
	local evap_max=r(max) //Get max value
	replace new_swp_uec=new_swp_uec/`evap_max' //Normalize UEC bt 0 and 1
betafit new_swp_uec //Fit beta distribution
	local evap_alp=e(alpha)
	local evap_bet=e(beta)
dpplot new_swp_uec, dist(beta) param(`e(alpha)' `e(beta)')

*Water Heater
sum new_wht_uec  
	local wht_max=r(max) //Get max value
	replace new_wht_uec=new_wht_uec/`wht_max' //Normalize UEC bt 0 and 1
betafit new_wht_uec //Fit beta distribution
	local wht_alp=e(alpha)
	local wht_bet=e(beta) 
	
*Clothes washer
sum new_cws_uec  
	local cws_max=r(max) //Get max value
	replace new_cws_uec=new_cws_uec/`cws_max' //Normalize UEC bt 0 and 1
betafit new_cws_uec //Fit beta distribution
	local cws_alp=e(alpha)
	local cws_bet=e(beta) 	
	
*Dryer
sum new_edy_uec  
	local edy_max=r(max) //Get max value
	replace new_edy_uec=new_edy_uec/`edy_max' //Normalize UEC bt 0 and 1
betafit new_edy_uec //Fit beta distribution
	local edy_alp=e(alpha)
	local edy_bet=e(beta) 	
	
*Dishwasher
sum new_dwh_uec  
	local dwh_max=r(max) //Get max value
	replace new_dwh_uec=new_dwh_uec/`dwh_max' //Normalize UEC bt 0 and 1
betafit new_dwh_uec //Fit beta distribution
	local dwh_alp=e(alpha)
	local dwh_bet=e(beta) 	
	
*Pool pump
sum new_pmp_uec  
	local pmp_max=r(max) //Get max value
	replace new_pmp_uec=new_pmp_uec/`pmp_max' //Normalize UEC bt 0 and 1
betafit new_pmp_uec //Fit beta distribution
	local pmp_alp=e(alpha)
	local pmp_bet=e(beta) 		
	
*Spa (other)
sum new_spa_uec  
	local spa_max=r(max) //Get max value
	replace new_spa_uec=new_spa_uec/`spa_max' //Normalize UEC bt 0 and 1
betafit new_spa_uec //Fit beta distribution
	local spa_alp=e(alpha)
	local spa_bet=e(beta) 		
	
*Spa (electric)
sum new_sph_uec  
	local sph_max=r(max) //Get max value
	replace new_sph_uec=new_sph_uec/`sph_max' //Normalize UEC bt 0 and 1
betafit new_sph_uec //Fit beta distribution
	local sph_alp=e(alpha)
	local sph_bet=e(beta) 		
	
*Well Pump
sum new_wpm_uec  
	local wpm_max=r(max) //Get max value
	replace new_wpm_uec=new_wpm_uec/`wpm_max' //Normalize UEC bt 0 and 1
betafit new_wpm_uec //Fit beta distribution
	local wpm_alp=e(alpha)
	local wpm_bet=e(beta) 		
	
	
********************************************************************
********************************************************************
********************************************************************
*SIMULATION MODEL
********************************************************************
********************************************************************
********************************************************************	
clear

forvalues i = 1(1)`sim'{
set obs 5046   //Number of households
gen hh=_n      //HH indicator
gen te=rnormal(0.029,0.0086) //Drawing treatment effect (Column 4 of Table 2, Panel B)

***************
*STEP 1: Drawing random sample of HHs, appliances, and UECs
* Notes: - For each appliance: 
*            (i)  Make N Bernoulli draws (0s and 1s) for appliance ownership
*                  based on estimated mean ownerships. 
*            (ii) Draw UECs based on estimated beta distribution. Assign 
*                  random UECs to each household/appliance. 

*Evaporative cooling
gen evap_cool=uniform()<=`evap_sat'
	gen evap_uec=rbeta(`evap_alp',`evap_bet')
	replace evap_uec=evap_uec*`evap_max' //Bringing back up to scale
	replace evap_cool=evap_cool*evap_uec 
	drop evap_uec
	
*Water Heater
gen wtr_heat=uniform()<=`wht_sat'
	gen wht_uec=rbeta(`wht_alp',`wht_bet')
	replace wht_uec=wht_uec*`wht_max' //Bringing back up to scale
	replace wtr_heat=wtr_heat*wht_uec 
	drop wht_uec

*Water Heater (Alternative - Every House Owns)
gen wtr_heat_alt=uniform()<=1
	gen wht_uec_alt=rbeta(`wht_alp',`wht_bet')
	replace wht_uec_alt=wht_uec_alt*`wht_max' //Bringing back up to scale
	replace wtr_heat_alt=wtr_heat_alt*wht_uec_alt 
	drop wht_uec_alt
	
*Clothes washer
gen clothes_wash=uniform()<=`cws_sat'
	gen cws_uec=rbeta(`cws_alp',`cws_bet')
	replace cws_uec=cws_uec*`cws_max' //Bringing back up to scale
	replace clothes_wash=clothes_wash*cws_uec 
	drop cws_uec

*Dryer
gen dryer=uniform()<`edy_sat'
	gen edy_uec=rbeta(`edy_alp',`edy_bet')
	replace edy_uec=edy_uec*`edy_max' //Bringing back up to scale
	replace dryer=dryer*edy_uec 
	drop edy_uec	
	 
*Dishwasher
gen dish_wash=uniform()<=`dwh_sat'
	gen dwh_uec=rbeta(`dwh_alp',`dwh_bet')
	replace dwh_uec=dwh_uec*`dwh_max' //Bringing back up to scale
	replace dish_wash=dish_wash*dwh_uec 
	drop dwh_uec	
	 
*Pool pump
gen pool_pump=uniform()<=`pmp_sat'
	gen pmp_uec=rbeta(`pmp_alp',`pmp_bet')
	replace pmp_uec=pmp_uec*`pmp_max' //Bringing back up to scale
	replace pool_pump=pool_pump*pmp_uec 
	drop pmp_uec	
	 
*Spa (other)
gen spa=uniform()<=`spa_sat'
	gen spa_uec=rbeta(`spa_alp',`spa_bet')
	replace spa_uec=spa_uec*`spa_max' //Bringing back up to scale
	replace spa=spa*spa_uec 
	drop spa_uec	

*Spa (electric)
gen sph=uniform()<=`sph_sat'
	gen sph_uec=rbeta(`sph_alp',`sph_bet')
	replace sph_uec=sph_uec*`sph_max' //Bringing back up to scale
	replace sph=sph*sph_uec 
	drop sph_uec
	
*Well Pump
gen well_pump=uniform()<=`wpm_sat'
	gen wpm_uec=rbeta(`wpm_alp',`wpm_bet')
	replace wpm_uec=wpm_uec*`wpm_max' //Bringing back up to scale
	replace well_pump=well_pump*wpm_uec 
	drop wpm_uec	

********
*Scenario 1: Reduction Across all Water/Energy Appliances 
*Notes: - Assumes proportionate effect, i.e., 3% reduction in water use leads
*         to 3% reduction in electricity use. 
egen save_s1=rowtotal(evap_cool wtr_heat clothes_wash dryer ///
					 dish_wash pool_pump spa sph well_pump)
	replace save_s1=save_s1*te

sum save_s1
	scalar s1 = r(mean)

********
*Scenario 2: Reductions in Clothes Washer, Dishwasher, and Dryer  Only
*Notes: - Assumes proportionate decrease in water heater and well pump use. 
*       - According to Home Water Works (http://www.home-water-works.org/)
*         an old CW uses 40-45 gals/load and a HE CW uses 14-25 gals/load. 
*         Assume 35 gals/load here. 
*       - According to Home Water Works (http://www.home-water-works.org/)
*         older DWs use 10-15 gals/load and newers ones use less than 5.5 
*         gals/load.  Assume 8 gals/load here. 
*       - According to RASS data, households run approximately 256 
*         (se=166.37) clothes washer loads per year, 228 (sd=130.06) dryer 
*         loads per year, and 133 (sd=111.2181) dishwasher loads per year.  
*       - Multiplying CW and DW use, we get 256*35+133*8=10,024 gals/year. 
*       - Estimated water TE is 3,925 gals/year, so if all came off DW and 
*         CW this would be a 39.2% reduction in use. Assume same for dryer.       
egen save_s2a=rowtotal(wtr_heat well_pump)	
	replace save_s2a=save_s2a*te
egen save_s2b=rowtotal(clothes_wash dryer dish_wash)	
	replace save_s2b=save_s2b*0.392
gen save_s2=save_s2a+save_s2b
	drop save_s2a save_s2b

sum save_s2
	scalar s2 = r(mean)
	
********
*Scenario 3: Reductions in Pool & Spa Use Only
*Notes: - Assumes every house that owns a pool and spa stopped filling them,
*         i.e., let them evaporate, and stopped using their pool pump and 
*         heaters (100% reduction in electricity use for these).   
*       - Assumes proportionate decrease in well pump use. 
egen save_s3a=rowtotal(well_pump)	
	replace save_s3a=save_s3a*te
egen save_s3b=rowtotal(pool_pump spa sph)
gen save_s3=save_s3a+save_s3b
	drop save_s3a save_s3b	

sum save_s3
	scalar s3 = r(mean)
	
********
*Scenario 4: Every HH has an electric hot water heater
*Notes: - Assumes every house that owns a pool and spa stopped filling them,
*         i.e., let them evaporate, and stopped using their pool pump and 
*         heaters (100% reduction in electricity use for these).   
*       - Assumes proportionate decrease in well pump use. 
egen save_s4=rowtotal(evap_cool wtr_heat_alt clothes_wash dryer ///
					 dish_wash pool_pump spa sph well_pump)
	replace save_s4=save_s4*te

sum save_s4
	scalar s4 = r(mean)

*Saving estimated electric savings under each scenario
gen s1=scalar(s1) in 1
gen s2=scalar(s2) in 1
gen s3=scalar(s3) in 1
gen s4=scalar(s4) in 1
	drop if missing(s1)
	
drop save_s1 save_s2 save_s3 save_s4 hh evap_cool wtr_heat ///
     clothes_wash dryer dish_wash pool_pump spa sph well_pump te wtr_heat_alt
	 
save $temp\save_`i'.dta, replace
drop s1	s2 s3 s4 
}
*
use $temp\save_1.dta, clear
forvalues i = 2(1)`sim' {
append using $temp\save_`i'.dta
}
* 
forvalues i = 1(1)`sim' {
erase $temp\save_`i'.dta
}
* 
*Labeling and savings results
 //Note: Estimated yearly TE is 297.84/year 
label var s1 "Scenario 1 - All Reduction"
label var s2 "Scenario 2 - CW, DW, and Dryer Reduction"
label var s3 "Scenario 4 - Pool Reduction"
label var s4 "Scenario 5 - All Reduction: 100% Electric Wtr Heater"
save $clean\sim_results_$outputdate, replace

frame change data_main
frame drop simulation